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ABSTRACT 

Remnant planetesimals might have played an important role in reducing the 
orbital eccentricities of the terrestrial planets after their formation via giant im- 
pacts. However, the population and the size distribution of remnant planetesi- 
mals during and after the giant impact stage are unknown, because simulations 
of planetary accretion in the runaway growth and giant impact stages have been 
conducted independently. Here we report results of direct N-body simulations of 
the formation of terrestrial planets beginning with a compact planetesimal disk. 
The initial planetesimal disk has a total mass and angular momentum as ob- 
served for the terrestrial planets, and we vary the width (0.3 and 0.5AU) and the 
number of planetesimals (1000-5000). This initial configuration generally gives 
rise to three final planets of similar size, and sometimes a fourth small planet 
forms near the location of Mars. Since a sufficient number of planetesimals re- 
mains, even after the giant impact phase, the final orbital eccentricities are as 
small as those of the Earth and Venus. 



Subject headings: Accretion, terrestrial planets 
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Introduction 



The rocky planets are believed to have formed via the accretion of small planetesimals. 
The formation mechanism of planetesimals and their initial mass and spatial distribution 
are however still controversial. The standard picture of accretion of terrestrial planets from 
planetesimals is as follows. During the early stages of planetesimal ac cretion, larger planetes- 



imals grow faster than sm aller ones owing to their stronger gravity (IGreenberg et al.lll978 



Wetherill fc Stewartl Il989l ). Through this runaway growt h stage, a few tens of Ma rs-size 
protoplanets form with mutual separation of 10 Hill radii (IKokubo fc Idalll998l . 120021 ). The 
growth time scale of protoplanets i s estimated to be ~ 0.1-lMyr around 1AU, and is longer 
with larger distance from the Sun (j Wetherill fc Stewartl 1 19931 ; IKokubo fc Idall2002l ). As long 
as some amount of nebular gas and/or planetesimals remain, its damping effect stabilizes 



the o rbits of protoplanets, preventing mutual collisions (jlwasaki et al.ll2002 



Kominami Ida 



20021 ). As the amount of remnant gas and/or planetesimals decreases, the orbital eccentrici- 
ties of protoplanets increases due to their mutual interactions. Eventually their orbits become 
chaotic and late time giant impacts occur. During this giant impact st age, whose time scale is 
considered to be ~ lOOMyr, the current terre strial planets form (e.g. IChambers fc Wetherill 



19981 ; lAgnor et al.lll999t IKokubo et al.ll2006l ). The orbital eccentricities of planets immedi- 
ately after giant impacts are likely to be much larger than those of the current terrestrial 
planets. Therefore, interactions with remnant gas and/or planetesimals is expe cted to re- 
duce t h eir eccentricities . Revi ews for the processes described above are given by IChambers 
tooi ): iNagasawa et al.l J2007h . 



Whether remnant gas or remnant planetesimals is more important for reducing eccen- 
tricities primarily depends on the time scale of gas dissipation. If the time scale is long, 
the velocity dispersion of planetesimals is suppressed by the gas drag. Hence, the gravi- 



tational focusing effect of protoplanets is enhanced, resulting in a fas 
planetesimals and a lower eccentricity distribution of the final planets ( 


: clean up of remnant 


Aenor & Ward 


2002 


Kominami & Ida 


2002, 


2004; 


Naeasawa et al. 


2005; 


Ogihara et al. 


2007 


). On the other hand. 



if the time scale of gas dissipation is short, planetesimals remain unaccreted by protoplanets 
for a longer perio d of time. In th i s case, remnant planetesimals may be req uired to reduce 



the eccentricities (jChambersll200ll ; lO'brien et al.ll2006l ; iRaymond et al.ll2006l ). In this paper 



we examine the latter scenario. In other words, we ignore the effects of the gas drag and the 
tidal interaction between a gas disk and protoplanets. The effects of gas will be investigated 
in a future study. 

There have been several at tempts to examine the effect of remnant planetesimals base d 
on direct A^-body simulations (j Chambers! 2001; O'brien et al.ll2006t IRaymond et al.ll2006l ). 
and with simulations using a hybrid-code (IKenyon fc Bromley! 120061 ). Direct A^-body simu- 
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lations usually adopt lunar to Mars size protoplanets surrounded by smaller planetesimals as 
initial conditions. Direct iV-body simulations suggest that the eccentricities of final planets 
are further reduced as the total mass of planetesimals increases. Even for the same total mass 
of planetesimals, the damping effect is stron ger with a larger number of smaller planetesimals 
(jO'brien et al.ll2006l ; iRaymond et al.l 120061 ) . However, the total mass and mass distribution 
of the remaining planetesimals are unknown, unless mass evolution in the runaway stage is 
followed. 



The hybrid code of iKenyon fc Bromley! (120061 ) ; iBromley fc Kenyonl (120061 ) . is able to 
follow planetary accretion through both the runaway and giant impact stages. In their code, 
the mass and velocity distributions of planetesimals contained in multi annuli are solved 
by a statistical approach whereas orbits of protoplanets are solved using direct iV-body 
calculations, that include the effect of interaction between planetesimals and protoplanets. 
Statistical approaches based on the local approximation produce cons istent results wit h 
those obtained from direct iV-body simulations in the runaway stage (e.g. llnaba et al.ll200ll ). 
However, it is questionable if statistical approaches can accurately follow the late stage of 
planetary accretion because orbital eccentricities of remnant planetesimals are usually very 
large. 

Here we report results of direct iV-body simulations beginning with a planetesimal disk 
until the end of planetary accretion in the terrestrial region. We consider compact planetesi- 
mal disks (initial disk widths of < 0.5AU), whose total masses and total angular momenta are 
the same as those of the present terrestrial planets. These initial conditions are used since ac- 
cretion simulations beginning from compact disks are usually computationally less expensive 
than those from wider disks. This is the case even with the same initial number of particles, 
because of the rapid decrease in the number of particles through accretion. Another reason 
to adopt compact disks is that total angular momenta of final planetary systems obtained 
from most of previous simulations are much larger than for the terrestrial planets, as thes e 



simulations usually have a super- massive Mars (e.g. I Chamb erj |2 00 It IRaymond et al.ll2006l ). 
This excess angular momentum is likely due to initially extended disks. Though Jupiter 
removes angular momentum, mostly from the asteroid region, its effect does not seem to 
be sufficiently strong in the terrestrial region. As one possibility for this issue, we consider 
initially compact disks, supposing that they result fro m, for example, dust migration due to 



the gas drag prior to formation of planetesimals (e.g. lYoudin fc Shull2002l ). 



In § 2, we explain the numerical methods used in this study. We show results of simula- 
tions in § 3. We compare our results with previous simulations in § 4. In § 5, we give some 
physical interpretations for our simulation results using analytic estimations. We summarize 
our results in 5 6. 
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Methods 



The runaway growth stage is shorter than the giant impact stage, but orbits of a large 
number of bodies need to be followed. On the other hand, although fewer bodies are necessary 
for the giant impact stage (unless the effect of fragmentation is considered) , more care must 
be taken to accurately follow the orbital evolution over many more dynamical times. Taking 
these physically different types of the accretion stages into account, we apply different N- 
body codes to these two stages of the evolution. 



The runaway g rowth stage is simulated with the parallel tree-code PKD CiMF ( iRichardson et al 
2000l ; IStadell l200ll ) for 10 5 yr with adopting artificially enhanced radii. The code uses a 



fourth-order multi-pole expansion for the force calculations, and a second-order leap-frog 
scheme is used for time integration. We apply a hierarchical time stepping with the largest 
time step of 1.8 days (0.005 yr). The opening angle of 0.5 is used as a criterion for search- 
ing down the tree. The energy error in the runaway stage is \ AE/E\ ~ 10~ 4 -10~ 3 , arising 
entirely from the integrator. The error due to the force calculation using our tree method is 
negligibly small. Using the output of the runaway stage as the i nitial condition, w e simulate 
the giant impact stage with the hybrid symplectic code Mercury ( Chambe"rslll999l ) for 2 x 10 8 
yr witho ut any enhancement of radii. This code uses a mixed variable symplec tic (MVS) 
method ( Kinoshita et al. 1991 ; Wisdom fc Holman 1991 ; Saha fc Tremainej|l992l ) for orbits 
around the Sun whereas close encounters are integrated by the Bulirsch-Stoer method. We 
use a fixed time step of 6 day s, which is the s ame as or similar to those adopted in the 
previous works using Mercury (jChambersI l200ll ; iRaymond et al.l 120061 ) . The energy error in 
the giant impact stage is \AE/E\ ~ 10~ 5 . It usually takes less than one computer day for 
a simulation of the runaway stage with PKDGRAV, whereas it can take several months to 
compute the giant impact stage with Mercury. 

Whereas MVS type integrators can take much larger time steps than those used with 
the leap-frog scheme, the Mercury code uses direct summation for calculations of the mutual 
gravity force. PKDGRAV is thus faster than Mercury as long as the number of particles is 
larger tha n seve ral hundred. Some comparisons between these two codes are also found in 
Raymond! J2005h . 



The enhancement of radii in the runaway stage is used in order to reduce computational 
time; we use an enhancement factor of radii g = 4.3. This gives an analogous effect of the 
gas drag, and the growth time scale of protoplanets is reduced by a factor of ~ g 2 (Kokubo & 
Ida 1996, 2002, see also eq. [20]). On the other hand, the growth time scale of protoplanets 
is actually reduced by a factor of ~ 2 by the gas d rag, where (3 is the factor for reduction 
of planetesimal eccentricities (IKokubo fc Idal 120001 ) . Thus, our simulations approximately 
mimic a situation in which the gas disappears suddenly at (g/ ' (3) 2 x 10 5 yr. This time scale is 
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proba bly shorter than the typical life time of circumsteller disks, ~ a few Myr (IHaisch et al. 
20011 ) . although the exact time scale for formation of planetesimals from dusty gaseous disks is 
not known. As long as the number of planetesimals is sufficient, the orbits of protoplanets are 
stabilized during the runaway stage b y dynamical friction such that the growth mode is not 
affected by the enhancement of radii (IKokubo fc Ida! 1 1 9961 ). However, in the transition from 
the runaway growth stage to the giant impact stage, faster clean up of remnant planetesimals, 
due to this approach, usually causes final planets to have higher eccentricities. Thus the radii 
of particles should be set to realistic values before planetesimals are too depleted in order to 
accurately examine the effect of remnant planetesimals. Additional simulations and analytic 
calculations were performed to assess the sensitivity of the results to varying the time at 
which g is reduced to unity. 

We use 10 different initial conditions which are summarized in Table 1. The total mass 
and angular momentum are assumed to be the same as for the present terrestrial planets 
(1.98M E and 1.86M E AU 1/2 y/UM^, respectively, where M E , G, and M stand for the mass 
of Earth, the gravitational constant, and the solar mass, respectively) with the central star's 
mass equal to the solar mass. The initial width of a planetesimals disk A^k is taken to 
be 0.3 AU or 0.5AU. The planetesimal mass is assumed to be identical and the number of 
planetesimals N varies from 1000 to 5000. The physical density of all the bodies is assumed 
to be p = 2g cm -3 . The surface number density n(a) as a function of the semimajor axis 
a is given by a power law n(a) oc a a with a = — 1 or —2. We also conduct two additional 
simulations for N = 1000, where we switch the code and g at 5 x 10 4 yr in order to check 
whether the outcomes are affected by this timing. 

It would be very interesting to investigate the accretion of planets using more extended 
disks (A disk > 0.5AU), but computationally too expensive with our current codes. The rate 
at which planetesimals merge is slower in a wider disk, particularly at its extremities, so we 
need to use the tree method for the gravity calculation for longer period of time. On the other 
hand, it is not appropriate to use the leap-frog integrator for the long-term orbital evolution 
for the following reasons. Firs tly, since the leap-frog i ntegrator causes a secular error in 
the longitude of the perihelion (IKokubo fe Makind 120041 ). it does not accurately treat long- 
term secular interactions. Secondly, since the standard (or explicit) block multi- timestep 
algorithm used for the leap-frog int egrator is not ti me- symmetric, the error in the energy 



accumulates with close en counters (IHut et al. 
order Hermite- integrator ( Kokubo fc Makinol 



1995). This is also the case for the higher 
20041 ). We ensure that the transition from 



PKDGRAV is chosen conservatively, thus we achieve high energy conservation as mentioned 
above. 



Although the implicit block time-step algorithm can avoid this problem (jMakino et al 
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20061 ). MVS integrators have considerable advantages for simulations of plan etary accretion. 



We are therefore implementing the SyMBA integrator (IDuncan et al.lll998l ) into the latest 
version of PKDGRAV, which enables us to simulate planetary accretion in wider disks. These 
simulation results will be reported later. 



3. Results 

3.1. An example of evolution: Run 6 

As an example, we first explain time evolution of Run 6. For this simulation the width of 
the initial disk Adisk is 0.5AU and the initial number of planetesimals N is 3000. Figure 1-3 
show time evolution of this simulation: snapshots on the plane of the semimajor axis versus 
the orbital eccentricity (Fig. 1), the cumulative number of planetesimals plotted against the 
mass (Fig. 2), and the epicyclic velocity plotted against the mass (Fig. 3). For detailed anal- 
ysis, we divide the accretional evolution into four different stages (the runaway, oligarchic, 
giant impact, and post giant impact stages), rather than two main stages discussed so far. 



3.1.1. Runaway growth stage (~ 10 4 yr) 



In the early stage, most of the mass of the system is contained in smallest planetesimals. 
In this case the epicyclic velocity, v = (aCly/ e 2 + i 2 ) (where Vt is the orbital frequency, 
and e and i are the orbital eccentricity and inclination of a planetesimal, respectively), is 
regulated by the smallest planetesimals and is typically as large as their escape velocity; 
w esCj o = ^/2Gm / (gr ) (Fig. 3), where m and r are the mass and the radius. If v is much 
smaller than the escape velocity of the largest body t> esc ,p = y/2Gm p / (gr p ) (where m p and 
r p are the mass and the radius respectively), and is a decreasing function with mass as 
shown in Figure 3, then the largest body s tarts to grow much faster than nearby objects . 



This growth mode is called runaway growth (jWetherill &: Stewartlll989l ; iKokubo &: Idalll996 
Weidenschilling et al.lll997l ). At t ~ 10 4 yr, the power-law index q (dN c oc m g dm, where A, 



ation 



1996 



is the cumulative number and m is the planetesimal mass) is about —2.7 in our simu 
(Fig. 2). This value is close to q ~ —2.5 obta ined in simula t ions o f Kokubo fc Ida! 
2000) and the analytical estimate q = -8/3 bv iMakino et all Jl998h . jMakino et all fll998f l 
assumed complete energy partitioning (v oc m -1 / 2 ) in the strong gravit ational limit ( v <C 



sec 



^esc,o)- However, the actual velocity distribution is less steep than this (jRafikovi 12003k 
also our Fig. 3). If we assume v oc m 7 with 7 ~ —1/4 , whic h is a rough approximation of 
Figure 3 at t ~ 10 4 yr, the formulation of IMakino et al.l (119981 ) gives q = —13/6 + 7 ~ —2.4.) 
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3.1.2. Oligarchic growth stage (~ 10 5 yrj 



As large bodies grow, their mutual gravitat ional interactions leave their orbits separated 
by 5-10 Hill radius (IKokubo & Idalll995l . Il998l . see also our Fig. 5). The Hill radius rn of a 
planetesimal of mass m p is given by 



r H = ahr 



2mp 

3M m 



1/3 



where a is its semimajor axis and h p is the reduced Hill radius. The largest bodies gravita- 
tionally influence the velocity evolution of all the neighboring planetesimals (llda &: Makino 
19931 ). increasing towards the escape velocity of the protoplanet v eS c,p ( see the panel of t — 10 5 
yr in Fig. 3). On the other hand, the velocities of protoplanets v p are also influenced by 
the energy partitioning with surrounding smaller planetesimals. Indeed, the value of v p for 
Run 6 is quite close to the equilibrium value v p>eq (~ v eBCj o), which is theoretically estimated 
neglecting the mutual perturbations of protoplanets (see eq. [35] in § 5.3). 

Since the growth rate of the largest body slows down at the expense of its nearby 
neighbors, the large st intermediate mass objects begin catch up with the largest bod y 
( llda fc Makinolll993l ). This growth mode is called oligarchic growth ( jKokubo fc Idalll998l ). 
At t — 10 5 yr in Run 6, about half of the total mass is contained in the 10 largest oligarchic 
bodies. Since the growth of smaller planetesimals has substantially stalled, protoplanets 
start to separate from the continuous size distribution. Therefore, the number of planetes- 
imals decreases mostly by accretion onto protoplanets and not by mutual collisions. Since 
v for planetesimals is nearly independent of the mass, so is their collision probability with 
protoplanets. Hence, the power-law index q for the mass distribution of planetesimals does 
not change from —2 after this stage. 



3.1.3. Giant impact stage (~ 10 6 yr) 



Without the damping force by remnant planetesimals and/or gas, a multiple proto- 
planet system undergoes an orbital instability after a certain time T inst . This instability 
tim e Tinst depends on the orbital separation, eccentricities, and absolute mass of protoplan- 



ets ( 


Chambers et al. 


1996: 


Ito & Tanikawa 


1999; 


Yoshinaea et al. 


1999; 


Iwasaki & Ohtsuki 


2006 


). We switch the code and reduce the value of g from 4.3 to 1 at t — 10 5 yr. For 



Run 6, T inst this time is estimated to be 10 5 -10 6 yr from above studies. With a decreasing 
total mass of planetesimals that have higher velocities, the damping due to the dynamical 
friction of planetesimals becomes less effective. At this point the orbital instability and mu- 
tual collisions of protoplanets start to take place. We find that the orbital instability starts 
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immediately after t = 10 5 yr, and the number of protoplanets decreases from 12 at t = 10 5 
yr to 10 at t = 2 x 10 5 yr (here we assume a protoplanet to be a body with mass > 2 x 10 26 g 
~ 50m ). 

During the giant impact stage, the mass distribution changes mainly owing to collisions 
between protoplanets, while the population of small planetesimals does not change so much. 
This can be seen in Figure 2; from t = 10 5 yr to t = 10 6 yr, when the number of protoplanets 
reduces from 12 to 5, whereas the total number of particles reduces only from ~ 400 to 
~ 300. Because the velocities of protoplanets are much smaller than those of planetesimals, 
mutual collisions of protoplanets occur quickly. This is similar to results of simulations 
including the damping f orce due to the tidal interaction between a gas disk and protoplanets; 



Kominami &: Idal (120021 ) find that the giant impact stage becomes shorter with the stronger 



damping force. 



3.1.4- Post giant impact stage (> 10 7 yr) 

The number of protoplanets further reduces to 3 at t — 10 7 yr, after two final giant 
impacts that occur shortly before t = 10 7 yr. Through the giant impact stage, the mutual 
separation between protoplanets normalized by their Hill radii increases to ~ 30. The mass 
distribution becomes completely bimodal (Fig. 2) with the masses of protoplanets smaller 
than the isolation mass by a factor of ~ 3. The isolation mass is the total mass contained in 
a ring of width 30 Hill radius with the initial surface density (see eq. |19|). This deviation 
likely comes from decrease of the surface density (by a factor of ~ 2) due to expansion of 
the disk from its initial diameter via gravitational scattering of protoplanets. For remnant 
planetesimals, the power-law index q remains to be ~ —2 and the largest mass is ~ 50m , 
which is similar to the protoplanet 's mass during the runaway to oligarchic stages. 

Since the mutual interactions between protoplanets after the giant impact stage is rather 
weak, the eccentricities of protoplanets are expected to be determined by the energy parti- 
tioning with remnant planetesimals. As we discussed for the oligarchic stage, v p in the giant 
impact phase is surprisingly close to the equilibrium value (~ t> eS c,o; see ec L- although 
some amount of remnant planetesimals may be necessary to achieve full equilibrium. For this 
simulation we find that the fraction of the total mass contained in planetesimals is 0.29 and 
0.16 at t — 10 6 and 10 7 yr, respectively. Naively, one might predict that equilibrium occurs 
once the mass in planetesimals is comparable to the total mass of protoplanets (see § 5.2). 
This slight contradiction might mean that the damping due to giant impacts themselves or 
some other unknown mechanism, works effectively. 
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3.2. Evolution of orbital spacing and eccentricities 

Here we analyze the evolution of each simulation quantitatively and discuss its depen- 
dence on the initial conditions. In order to examine characteristics of the evolution of the 
largest bodies, we define a planet (or a protoplanet) as a body having mass larger than 
2 x 10 26 g. The physical meaning of this choice is that a body larger than this mass regulates 
the velocity evolution of all the neighboring bodies (see § 5.3). In this case the following 
discussion does not strongly depends on the definition of the minimum mass of planets. We 
will explore the following four quantities: [1] the number of planets, N p , [2] the mass fraction 
of planets compared to the total mass, f p , [3] the orbital spacing of planets normalized by 
the mutual Hill radius, b p , [4] the eccentricity of planets, e p , and [5] the ratio of the effective 
mass of planetesimals to the mean mass of protoplanets m e ff/(m p ). The effective mass of 
planetesimals is defined as m e g = (m 2 ) / (m), where (m 2 ) and (m) are the mean squared and 
mean masses of planetesimals, whose masses are smaller than 2 x 10 26 g. The mass ratio 
^eff/(^ P ) would be important for the evolution of e p , as its equilibrium value due to the 
dynamical friction is given by (eq. [T5] in § 5.1) 



p,eq — 



4m cff 
3(m p ) 



(e 2 ) 



2\l/2 



(2) 



We apply the following form for the averaged orbital spacing normalized by the mutual 
Hill radius 



with the reduced mass 



and the mutual hill radius 



\ 



EjV(( Q i+i - Q j)A"Hj) 2 /^ 



i'j 



m p,j m p,j+i 
m Pj + m Pj+i : 



(3) 



(4) 



r H,. 



{dj + a j+1 ) 



m p,j + m P,j+i 
3M. 



1/3 



(5) 



Here aj , ej , and Tflp j £1X6 the semimajor axis, orbital eccentricity, and mass of the protoplanet 
j in the order of semimajor axis, respectively. We use the following form for the averaged 
eccentricity, which characterizes the energy of epicyclic motion of planets: 



EiVp 2 
j=l m p,j e j 



(6) 
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where e,- is the orbital eccentricity of the protop lanet j. The mass weighted eccentricity (e.g 



Bromley &: Kenyonll2006l ; iRaymond et al.ll2006l ) gives a slightly smaller value than the above 



eccentricity. 

Figure 4 shows the evolution of these four quantities for the case of Adisk = 0.3AU. 
In the third panel from the top, we also plot the orbital instability time as a function of 



the or bital spacing for the two cases of e p /h p = 2 and 4, respectively, from lYoshinaga et al. 



(119991 ) (see eq. 23). Here h p = (2(m p )/(3M )) 1 / 3 is the averaged reduced hill radius. These 
lines indicate the stability of multiple protoplanet systems; if the orbital spacing is narrower 
than these lines, orbital instability will occur. Several proto planets form at y 10 4 yr and the 



normalized orbital spacing b p is about 10, as pointed out by lKokubo fc Idal (119981 ). Since the 
orbital instability time is 10 5 -10 6 yr in such systems, giant impacts between protoplanets 
start around that time as we discussed in the previous section. The orbital eccentricity prior 
to the giant impact stage is ~ 0.03 while it increases up to 0.1 during the giant impact stage. 
The corresponding normalized eccentricities, e p /h p , are 2-3 and ~ 10, respectively. Through 
the giant impact stage, b p increases to 20-30. As the radial excursion of planets during the 
giant impact stage determines the final separation of planets, the relation between b p in 
the final state and e p d uring the giant impact stage can be approximately represented by 



(IKominami fc Idall2002h 



b p h p ~ 2e p . (7) 

This is roughly consistent with our simulation results. 

During the post giant impact stage, the orbital eccentricities are reduced (e p /h p ~ 3) and 
these values have little dependence on initial parameters (we discuss the weak dependence 
in detail in § 5.3). Since b p is large enough in the final state, the mutual interaction of 
planets is likely to be unimportant. In this case, final eccentricities of planets are expected 
to be determined by the energy partitioning with remnant planetesimals. The mass ratio 
^eff/(^p) decreases nearly monotonically with time and is about 0.01-0.02 in the final 
state (the bottom panel of Fig. 4; the decrease of m e s/{m p ) is due to the increase of (m p ), 
while m e g is nearly constant except at very early times). Since the mean eccentricity of 
planetesimals is ~ 0.3 in the post giant impact stage (Fig. 3), the equilibrium eccentricity of 
planets (eq. [2]) is estimated to be 0.03-0.05. This is almost the same as the values obtained 
in our simulations. 

While the final giant impact occurs before t ~ 10 7 yr in most of the runs, it happens 
at t ~ 5 x 10 7 yr for Run 1 (N = 1000). Because of that impact, the orbital spacing 
for Run 1 becomes even wider. This unstable behavior likely suggests that the dynamical 
friction for N = 1000 is less effective as compared with larger N. However, except for this 
event, the dependence of the evolutions of all the quantities shown in Figure 4 on N is very 
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small. Switching g in the earlier time does not affect the quantities in the final state shown 
in Figure 4, either, although Run lb is dynamically more excited during the giant impact 
stage. This is because the mass variation of protoplanets in Run lb is somewhat large and 
smaller protoplanets are dynamically enhanced by larger protoplanets. 

Figure 5 is the same as Figure 4, but for the case of Adisk = 0.5AU. The evolution 
of all the quantities are very similar to those in Figure 4, although the early evolution is 
slightly slower due to the lower initial surface density. The dependence on N is very small 
here as well as in Figure 4. In particular, the final b p has converged to ~ 30 for all the 
runs. However, the final eccentricity e p is much larger for the case with earlier switching 
of g (Run 5b). This is in an opposite sense to our prediction because dynamical friction 
works more effectively with a higher density of planetesimals. We will discuss this issue in 
the next section along with the final configurations of the systems. Except for Run 5b, e p 
and m e ff/(m p ) in Figure 5 are slightly smaller than those in Figure 4. The smaller m e fj is 
due to the smaller surface density (the relation is roughly given by m e g- oc X 3 / 2 , where £ is 
the initial surface density of planetesimals; see § 5.3). 



3.3. Final systems 



Here we present the orbital parameters of all the final systems obtained in our simula- 
tions. Figure 6 shows snapshots of all of our runs on the a-e plane at 200Myr. Also, the 
number of planets N p , the averaged eccentricities of planets e p , and the angular momen- 



o tJH at 

(Laskar 


1997; 


iyi die ouu 

Chambers 


2001) 



s d 



i m p,.r 



1 



1 - # 



COS Ij 



where ij is the orbital inclination of the planet j. We take 5Myr averages for e p and Sd- 
For the current terrestrial planets, we take the mean values between the minim um and max- 



i mum orbital eccentricities and inclinations from 3Myr orbital integrations in iQuinn et al. 



( 11991I ). This operation roughly corresponds to subtracting only the free eccentricities, pro- 
vide d that the free eccentricit ies is larger than the forced eccentricities due to giant planets 



sec 



Murray fc Dermott 



tions obtained in lQuinn et al. 



1999 



, Chap . 7.4). In fact, the minimum eccentricities and inclina- 
(I199ll ) are almost zero except for Mercury. This suggests that 
the free eccentricity (inclination) and the forced eccentricity (inclination) are comparable for 
the current terrestrial planets. 



We always obtain three similar size planets between 0.5AU and 1.3 AU, except Run 1 has 
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only two planets in this region. The orbital spacing between planets are also quite similar. 
We find that more than 95% of the initial mass and angular momentum are contained in 
planets at 200Myr in all our simulations. Both the averaged eccentricities and the angular 
momentum deficits obtained from our simulations are comparable or even smaller than those 
for the current solar system, except for Run 5b and Run 8 (Table 1). Except these two runs, 
both e p and Sd are smaller for Adisk = 0.5AU than for Adisk = 0.3AU. There seems to be 
a weak tendency that e p and Sd decrease with increasing N. These trends are interpreted 
in terms of the effective planetesimal mass m eS , if the final e p is determined by the energy 
partitioning with planetesimals neglecting the mutual interaction of planets (§ 5.3). 

Differing from other runs, the mutual interaction between the innermost two planets in 
Run 5b and Run 8 is important even at the end of simulations, as their orbital separation is 
narrow. This seems related to the spatial distribution of planetesimals during the giant im- 
pact stage. In these runs, we find that two innermost planets quickly sweep out planetesimals 
in the inner region whilst there are still large numbers of planetesimals in the outer region. 
The outer planets with these planetesimals tend to push the middle planet inward. As a 
result, the two innermost planets continue interacting without sufficient dynamical friction 
due to surrounding planetesimals. Although our statistics is not sufficient, such a difference 
in the inner and outer region would tend to appear when g is reduced at earlier times or if 
the initial planetesimal mass were concentrated in the inner region. 

This fact seems related to the existence of the small outermost planet. When we compare 
simulations Run 1 and Run lb, the sizes and locations of two largest planets are very similar. 
While Run lb has inner and outer small planets, Run 1 has only an outer small planet (near 
2.1AU). Similarly, while the sizes and locations of the three largest planets in Run 5 and 
Run 5b are similar, only Run 5b has an outermost planet. These facts suggest that the 
tendency to have small planets in the inner and outer edges is stronger in simulations with 
earlier switching of g. We interpret this as due to a larger amount of planetesimals that are 
scattered inward (outward) at the inner (outer) edge of the disk before they are accreted 
by planets. Small planets form from these scattered planetesimals. A similar trend is also 
found in the simulations starting with a stronger gradient of the surface density in the radial 
direction (Runs 4 and 8). In these systems, the inner planets form quickly while large planets 
have not grown in the outer region. Then inner planets gravitationally scatter planetesimals 
and small protoplanets outward. The orbital eccentricities of protoplanets scattered outward 
are reduced by the dynamical friction of similarly scattered planetesimals. Eventually, these 
protoplanets can have stable orbits near the location of Mars and slightly grow as they collide 
with planetesimals. 
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4. Comparison with previous works 



Our simulation results are different from those starting with only protoplanets. iKokubo et al 

(120061 ) conducted simulations starting with ~ 15 mars-size protoplanets at 0.5 < a < 1.5AU 
and found that most of the final mass is contained in the largest two planets. Their fi- 
nal orbital eccentricities are usually higher than those for t he current terr e strial planets. 
Chambers fc Wetherilll (119981 ) also showed similar results to IKokubo et al.l (120061 ) for the 
case without perturbation of Jovian planets (their Model A). They also found that mass 
concentration within fewer planets is strengthened by the presence of Jovian planets (their 
Model B). The difference between their results and ours suggests that the number of final 
planets increases with a stronger damping force, which makes the radial excursion of pro- 
toplanets and thus the final s eparation between neigh boring planets narrower. Indeed, the 
same trend was also found by lKominami fc Ida! (120021 ). who examined the effect of damping 
due to the gas disk on the accretion of protoplanets. 

In recent direct A-body simulations stating with planetary embryos with small planetes- 
imals, the total mass of planetesimals is half or less than hal f, and the orbita l separation of 



embryos is equal to or less than 10 in units of the Hill radii (I Chambers! l200ll ; lO'brien et al. 



20061 ; iRaymond et al.l 120061 ) . Since the dynamical friction of surrounding planetesimals is 
not strong enough to suppress the orbital instability with these initial conditions (see § 5.2), 
giant impacts start immediately before planets substantially grow by accretion. They adopt 
a nearly identical size distribution of planetesimals, which in principle, does not change in 
such enhanced systems. Therefore, the effective mass of planetesimals after giant impacts 
occur is still given by the initial planetesimal mass. Thus, the dependence of the equi- 
librium eccentricity of planets on the initial planetesimal mass (e Pieq ~ y/ m / (m p ) (e 2 ) 1 ^ 2 ) 
is much stronger than we find. In some of their simulations, mutual interaction amongst 
the final planets is insigni ficant and the final e ccentricities seem to be close to the equilib- 
rium value. For example, lO'brien et al.l (120061 ) adopt the initial mass of planetesimals to 
be m ~ 1/400Me. In the late stages of their EJS (eccentric orbits of Jupiter and Saturn) 
simulations, most of the mass supplied to terrestrial planets is from small planetesimals with 
very high speed. If we convert their impact speed ~ 20km s _1 in the late stage to the ec- 
centricity around 1AU, it gives (e 2 ) 1 / 2 ~ 0.5. Supposing that the mean mass of planets is as 
large as the Earth's mass, we obtain e Pjeq ~ 0.025, which is even smaller than those for the 
current terrestrial planets, and consistent with their results. Therefore, we predict that the 
final eccentricities of planets would be further reduced if they adopted smaller planetesimal 
masses in their simulations. However, such a small effective mass of planetesimals in the 
beginning of the giant impact stage might be unlikely if we take the growth of planetesimals 
in the runaway and oligarchic stages into account. 



-14- 



Kenyon fc Bromleyl (120061 ) and iBromley fc Kenyonl (120061 ) conducted planetary accre- 
tion simulations starting with very small planetesimals (V p = l-5km) , using t heir h ybrid 
code. In their runs, simulations starting at 0.86-1.16AU in lKenyon fc Bromleyl (120061 ) have 
similar initial conditions to ours, although our disks are slightly more massive. The evolution 
of the number of oligarchic bodies (with masses >~ 10 25 -10 26 g in their simulations) and their 
orbital separation (their Hill parameter almost corresponds to 1/& P ) are very similar to our 
results. However, the final eccentricities of planets is more excited in their simulations (one 
of their simulations obtained three planets with e p ~ 0.1). In fact, in all of their simulations 
starting with wider initial disks (0.4-2.0AU), th e final planetary orbits ar e more eccentric 
than the current terrestrial planets (see Table I of lKenyon fc Bromleyll2006l ). and apparently 
remnant small planetesimals do not contribute to damping of eccentricities of planets. Since 
we have not conducted simulations starting with wide disks, it is not clear for us if their 
results obtained with a hybrid code are consistent with those obtained from direct A-body 
simulations. We are planning to conduct direct iV-body simulations with initially wide disks 
to clarify this problem. 



5. Analytic estimates 



In this section, we interpret our simulation results using analytic estimates. 



5.1. Evolution of velocities of planetesimals and protoplanets 



First we provide analytic formulation for the evolution of velocities of planetesimals 
and protoplanets necessary for subsequent discussions. Consider a situation in which proto- 
planets are spatially separated but mutually interact due to the distant perturbations and 
each protoplanet is surrounded by a swarm of planetesimals. We first consider equal-mass 
planetesimals, then the formulation is extended to the case of continuous size distribution 
(protoplanets are always assumed to be equal- mass). The mass, mean square eccentricity, 
and surface number density of planetesimals are represented as m, (e 2 ), and n, respectively. 
Corresponding characters for protoplanets are m p , (e 2 ), and n p , respectively. The scatter- 



ing cross section for planetesimal-planetesima 
planetesimal encounters m are given 



encounters a\ 



by Jlda fc Makino!ll993h 



and that for protoplanet- 



C P 



G(m + m') 



;afi) 2 ((e 2 ) + <e' 2 » 



(9) 
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m P -m' _ r ( g(m p + m') \ 2 

" \(am(el) + <e' 2 » J ' [W) 
where C e is the numerical factor of ~ 40, G is the gravitational constant, a is the distance 
of the system from the Sun, and Q is the orbital frequency. In equation ([91), we used primed 
characters, m' and (e /2 ), to distinguish two interacting planetesimal groups (the primed 
values are averaged later). Using the scattering cross sections, the change rates of (e 2 ) and 
(e 2 ) are given by 

1 ( e ) ( "' "\ _m-m'//„2\ , /„/2\\ , .. ( '"l> i ,>.,,-,.., / 2\ 



h§ - (d^) 2(4 " mB<e ' 2> - Wm ' <e ' )+ " mffl<e ' >)ffS " + ^' (12) 



where Ti nst is the time scale for orbital instability of protoplanets (see § 5.2). In equation (fTTj) . 
the first term stands for the viscous stirring due to planetesimal-planetesimal encounters 
while the second ter m for the viscous stirring due to protoplanet-planetesimal encounters 



( llda fc Makinolll993l ). In equation (fT2l) . the first term stands for the s um of viscous stirring 



and the dynamical friction both due to encounters with planetesimals (jlda fc Makinolll992l ) 
while the second term approximately accounts for the enhancement due to distant inter- 
actions between protoplanets. The rate of change of the inclinations are given by similar 
equations, but here we omit them. 

The surface number density per unit mass is given by dn/dm. As in the simplest case, 
we assume that (e /2 ) is independent of mass. In this case, after integration of equations (fTTl) 
and (fT2"|) over the range of m', nm' and nm' 2 in these equations can be replaced by 

m 'dn,= K, \m"in=^m a , (13) 

respectively. Here S s and m eff = (m /2 ) / (m f ) are the surface density and effective mass of 
planetesimals, respectively. Using these averaged quantities, we have simple implications 
from equations (TTT|) and (1T21) . Equation (TTTT) suggests that the velocity evolution of plan- 
etesimals is regulated by protoplanet-planetesimal encounters rather than by planetesimal- 
planetesimal encounters if 

n p m 2 p > / e S s m cfr , (14) 

where f e < 1 is the numerical factor associated with the velocity distribution. Equation (TL21 
suggests that the eccentricity of planets is given by 



3m p - m eff 
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in the equilibrium state provided that the mutual interaction between protoplanets is negli- 
gible. 

During the giant impact stage, in which orbits of protoplanets cross each other, the 
second term in the right hand side of equation (fl2|) can be replaced by the mutual viscous 
stirring term: 

n dt ) mnt 2 a ~ w [lb) 

with the cross section for the mutual scattering 

<®T* = C e ( j^hr V . (17) 



Kernel). 

The averaged eccentricity of protoplanets during the giant impact stage is determined by the 
balance between the mutual viscous stirring and the dynamical friction due to surrounding 
planetesimals. Assuming m p 3> m e s and m p (e 2 ) 3> m e e(e 2 ), we obtain 

<4>^=Q^J (e 2 ) 1/2 - (18) 



The same expression is obtained in iGoldreich et al.1 (120041 ) except for a factor of order unity. 



Equation (fl8l) somewhat overestimates (e 2 ) 1 ^ 2 during the giant impact stage as compared 
with those obtained from our simulations. This is probably because the motion of all the 
planets are not enhanced simultaneously in our simulations as some planets are in stable 
orbits isolated from others. 



5.2. Comparison of time scales and timing of the onset of giant impacts 

Here we discuss how the timing of the onset of the giant impact stage is affected by the 
radius enhancement factor g. Since the time and distance can be rescaled for our A-body 
simulations, g is physically associated with the ratio of the physical radius r to the Hill 
radius rn as g oc r/rn oc a _1 p -1 / 3 . 

For simplicity in this section, we call the most massive body in its feeding zone of 
width of & p rn, a planet (note that the definition of planets used in the main text follows the 
discussion in § 5.3). Defining the mass ratio of the planet with the total mass in the feeding 
zone to be f p , the planet mass m p is given as 

m p = 27r/ p a6 p r H S = (27r/ p 6 p S) 3/2 a 3 (2/3M ) 1 / 2 

Vioy V 2 0gcmy Vo.25y viau/ v ; 
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where £ is the surface density of all the bodies (thus assumed to b e constant regardless of 
/ p ). For f p = 1, m p corresponds to the so-called isolation mass (e.g. iKokubo &: Idall2000l ). 



Considering a two component system composed of planets and surrounding planetesi- 
mals, we introduce the following five timescales, which characterize the evolution of planets 
and surrounding planetesimals. [1] The growth time scale of planet T grow is given by 



T. 



grow 



1 dm r 



m p dt 



/ p )£cr co ifi' 



(20) 



with the collisional cross section ct co i (IGreenzewig fc Lissauerlll992l ): 



Ocol 



c C oi(gr p ) 2 



V 2 



(21) 



where c co \ is a numerical factor of ~ 8 (we assume (e 2 ) 1 / 2 = 2(z 2 ) 1 / 2 ) and v is the averaged 
velocity of planetesimals normalized by the escape velocity of the planet. [2] The depletion 
time scale of planetesimals Td ep is given by 



T, 



dep 



d(n p m p ) 



(1 - / p )E dt 



/ p ScT co if2 



(22) 



where n p = f p Y,/m p is the surface number density of planets. [3] The time scale for the 



evolution of the velocity of planetesimals T sca 
given by 

1 d(e 2 } 



due to gravitational scattering by planets is 



(e 2 ) dt 



/p^Osca^ 



(23) 



Here the scattering cross section cr sca corresponds to a s 
and (e 2 ) < (e /2 ) 



0" s 



c sca(fi ,? "p) ~ 4 



(eq. [TO] ) in the limit of m p ^> ml 



(24) 



where c sca is a numerical factor of ~ 16. Note that cr sca does not directly depend on g. 
However, as the velocity can be scaled by the escape velocity of the planet, which depends 
on g, a sca indirectly depends on g as represented by equation fl24l) . [4] The time scale for 
damping the velocity of planets Td amp due to dynamical friction of surrounding planetesimals 
is given by 

1 d{e 2 p ) ~ ! m r 



- damp 



(el) dt 



/p)S<j sca f2 



(25) 



[5] The time scale for the orbital instability T- mst of a multiple protoplanetary system with- 
out any damping force represents either the time of the first collision or the first close 
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nitions). The form of T inst is empirically given by 1 


Chambers et al. 


1996 


1999; 


Ito & Tanikawa 


1999: 


Iwasaki & Ohtsuki 


2006 


) 



log 



t, 



orb,l 



cib p + c 2 , 



(26) 



where T or b,i is the orbital period of the innermost protoplanet, and c\ and c 2 are numerical 
coefficients. These coefficients depend strongly on the orbital eccentricity and relatively 
weakly on the absolute averaged mass and the variation of masses. The depe n dence of c\ 
and c 2 on orbital eccentricities is summarized in Table III of lYoshinaga et al.l (119991 ). For 
b r 



10 and e p ~ 



4 fop, / iii^i 



t /T orbil ~ 10 5 - 10 6 . 
For simplicity, we normalize all the timescales as follows: 



T 



c co i£(gr p ) 2 



QT. 



(27) 



With this normalization, all the time scales except for Ti nst depend only on v and f p as 



f-i 

grow 



:w P )(~ + ^), 



^dep - /p( 3 +~2/' 



f-1 
damp 



2 (!-/p)^ 



(28) 



It should be noted that now the dependence on g is included only in Tinst- 



In order to obtain the time scales as functions of f p , we consider the evolution of v. 
In the early stages the smallest planetesimals dominate the mass of the system so that the 
velocity is as large as the escape velocity of planetesimals. The exact value of v at the initial 
state (when m p = m ) is determined by the balance between the mutual scattering and 
the collisional damping. Since the ratio of these time scales is given by the ratio of a co \ to 
cx sca , we obtain v = 1.17 for the initial state. Since the dimensional velocity of planets is 
1.17f esC)0 , v decreases as planets grow (or with increasing f p ). When f p reaches to a certain 
value, planets start to regulate the velocities of surrounding pla netesimals. In this c ase, v 
evolves as planets grow. Thus, v is obtained from T, 



grow 



T sca ( IDaisaka et al.l 120061 ) . The 



value of v increases with f p to 1.17 at f p = 0.5. For f p > 0.5, T^ ep becomes shorter than 
T sca , if v > 1.17. This means that planetesimals collide with planets before their velocities 
are further enhanced. Therefore, v takes a constant value, 1.17, for f p > 0.5. 
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We plot the evolution of v as a function of f p in Figure 7. In the same figure, we also 
plot T grow , T sca , T d amp, and T inst with the parameters used in Run 6 and b p = 10. For T inst , 
we take T inst /T or b = 5 x 10 5 with g — 1 and 4.3 as examples. This time scale would be 
appropriate judging from our simulations. It should be noted that T dep = T sca for f p > 0.5 
and T^ep further increases with decreasing f p for f p < 0.5, although we do not plot Td cp on 
Figure 7 in order to avoid confusion due to too many lines. 



Iwasaki et all (120021 ) and llwasaki fc Ohtsukil (120061 ) showed that orbital instability is 
preve nted when T darnp < CoT^t , with a coefficient c 3 of the order of unity (we assume c 3 = 3 
after llwasaki fc Ohtsukil (120061 )). Hence, the value of f p at the onset of giant impacts is 
estimated from the equation Td amp = 3T inst , and is 0.3-0.4 and 0.8-0.9 for g — 1 and 4.3, 
respectively, in Figure 7. However, Figure 7 also suggests that T inst (g = 4.3) > T dcp , which 
means that remnant planetesimals are depleted before the giant impact phase starts, as 
long as we keep g = 4.3. Complete dep l etion of pl anetesimals keeping b p ~ 10 is found in 
simulations with g = 6 in lKokubo fc Idal (120021 ) and lLeinhardt &: Richardson! (120051 ). and we 



also found the same results with additional tests. Therefore, in order to keep some amount 
of planetesimals at the onset of the giant impact stage, it must start before f p reaches 0.5 
(strictly speaking, this value is slightly higher than 0.5 for c 3 > 1). This also corresponds to 
the condition for the damping of enhanced eccentricities of planets during the giant impact 
stage. The condition under which the eccentricities of planets are damped before depletion 
of planetesimals is given by T da mp < Td ep , which gives f p < 0.5 (Fig. 7). For the case of 
g — 1, this condition is satisfied. The threshold value of g to satisfy the condition is roughly 
estimated to be 3. 

For our simulations (except Runs lb and 5b), we reduced g from 4.3 to 1 when f p is 
about 0.6 (Figs. 4 and 5). After reducing g, T inst becomes smaller than any other timescale 
(Fig. 7). Thus, the giant impact phase rapidly begins (the instantaneous reduction of g 
causes drop of v to ~ 1.17/\/4.3, but T sca is shorter than T grow by a factor of 4.3 with this 
small v and thus v increases near to 1.17 again before f p substantially increases). On the 
other hand, for Runs lb and 5b, f p is about 0.4 when we reduce g at 5 x 10 4 yr. In these cases 
giant impacts do not start immediately as Td am p ~ 3T inst at f p ~ 0.4 (this can be seen in 
Figs. 4 and 5 as N p does not change immediately after reducing g). In fact, after f p increases 
to ~ 0.5, giant impacts start in these simulations. Therefore, it is expected that reducing 
g at an earlier time does not affect the results (we have now conducted some simulations 
keeping fixed g = 1 and obtained consistent results; these results will be reported in another 
paper). 

To summarize, for simulations with a constant g throughout the entire accretion stage, 
sweeping up all the planetesimals and a subsequent giant impact stage are expected for 



-20- 



g > 3, whereas giant impacts during accretion of planetesimals and a subsequent damping 
of eccentricities of planets due to the dynamical friction are expected for g < 3. For g < 3, 
the evolution and final state of a system are expected to weakly depend on g, from the 
comparison between Runs 1 and lb and that between Runs 5 and 5b. 



5.3. Minimum oligarchic mass and equilibrium eccentricity of protoplanets 



Once the mass of the largest body is above a certain critical mass, the velocity evolution 
of neighboring planetesimals is primarily regulated by the largest body. We call this critical 
mass the minimum oligarchic mass. Then, the growth of smaller planetesimals near the 
largest body is stalled and as the largest body grows it starts to separate from the continuous 
size distribution of planetesimals. Therefore, the planetesimal size distribution after this 
stage is expected to be a continuous distribution with the maximum mass being the minimum 
oligarchic mass. This size distribution determines the equilibrium eccentricity of planets after 
the giant impact stage. Here we estimate the planetesimal size distribution, when the largest 
body starts to regulate the velocity evolution of all of its neighboring planetesimals. 

We consider the power-law size distribution for planetesimals dn oc m q dm with the 
upper and lower cutoff masses m p and m . The condition that the largest body regulates 
the velocity evolution is again given by (eq. 



ml > f e m T m eS , (29) 

where f e is a factor associated with the velocity distribution and is slightly smaller than 
unity and (= £ s /n p ) is the total mass of planetesimals (excluding the lar gest body or a 



plane t) in the heated region, where velocities are regulated by a planet (see llda fe Makino 
19931 ). with n p being the surface number density of planets (the sizes of the heated region 
and the feeding zones are similar). We also define the cumulative number in the heated 
region to be = n/n p , and now the size distribution is given by dN^ = km q dm. Using the 
condition for the mass of the largest planetesimals (or the second largest body in the heated 
region) rrih 

k q+1 

dN h = 2 mL ; , (for q < -1) (30) 

(31) 
(32) 



<m L -q-V 
the total mass and the effective mass m e g are, respectively, given by 

-q-2 



m T m eff 



"'L 

mo 



mdN h 



m 2 dN h 



, g + l 

'q + 2 

q + l 
g + 3 



1 - 



mo 



<?+3 



rrij 



(for q ^ -2) 



(for q 7^ —3) 
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For q = —2, we have mx = 2t«l \D.(mi,/mo). 

Substituting equation (j32J) and the relation m-L = 2 1// ^ +1 )m p into equation (j2"5j) . the 
condition approximately becomes 

1 < _ /e i±l 2 ^+ 3 )/^ +1 ), (33) 

g + 3 ' V ; 

which gives q > —2.2 for f e = 1. Since / e is expected to be slightly smaller than unity, q 
would be slightly smaller than —2.2 when the largest body starts t o regulate the v elocit y 



evolution. This might correspond to the value q ~ —2.5 obtained by lKokubo fe Ida! (119961 ). 
The fraction of planetary mass f p is given by 

^r^r^ m * /m ° ) - , - 2 - 1] ~ 1 - (34) 

J mo 11 ^ 

Substituting it into equation ffT9l) . we obtain the minimum oligarchic mass m P)0 ii = m p (q ~ 
—2.2) as function of q. This is shown in Figure 8 as well as the corresponding m e s for the 
parameters of Run 6. Using m p (q), we also plot the evolution of q as a function of f p in 
Figure 7. As we estimated above, Figure 7 suggests that v starts to increase at q ~ —2.5 
because of the gravitational scattering of planets. Therefore, when q ~ —2.5, growth of small 
planetesimals slows down as well as the evolution of their size distribution. 

However, at this stage, the velocities are still not high enough to suppress the growth of 
large planetesimals. Therfore, mass transfer from small planetesimals to large planetesimals 
further increases q to ~ —2. When q ~ —2, the velocity of planetesimals in the heated 
region is almost regulated by a single planet, whereas the contribution from a planet and 
all the other bodies are the same for q ~ —2.5. As a result, the velocity of planetesimals 
is enhanced toward the escape velocity of the planet. Then, actual separation of the planet 
from the continuous size distribution starts as we showed in § 3.1 (that is also why we keep 
q = — 2 for large f p in Fig. 7, although there is no physical reason for q to be strictly —2). 
Therefore, an appropriate minimum oligarchic mass to give the size distribution after the 
giant impact stage seems to be m p oli for q ~ — 2 (m po ii ~ 60m for q = — 2 whereas ~ 20m 
for q = —2.5). In the case of q = —2, m P)0 ii and m e s depend on the initial planetesimal mass 
m very weakly as m p>0 n(q = —2) oc [In (m P]0 ii/mo)]~ 3 ^ 2 S 3 / 2 and m e s(q = —2) = 0.5m PjO ii(g = 
-2)/ln(0.5m PiO ii/m ). 



Using the form of m e g for q = —2, we rewire the equilibrium eccentricity (eq. [2]) as 

,31n(0.5 mp , oB / mo )J UJ 2 >- (35) 



e P,eq 
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Substituting m Vi0 \i{q = —2), which is obtained from eqs. [19] and into eq. [35], we 
calculate e Pi6q for various cases of iV (or m ) and A disk (or £), assuming three final planets 
(thus m p = 0.66Me) and the velocity of planetesimals to be the escape velocity of a planet 

(v = v P . 



The calculated values are plotted in Fig. 9 and compared with our iV-body simulation 
results. We find good agreements between analytic estimates and simulations, as long as 
the planet-planet interactions are not important in the final state. Since the dependence of 
e Pi e q on N is very weak for N > 1000, we need more runs for statistics. We also conducted 
additional simulations for the case of iV = 100 and results are plotted in the same figure. 
In most of simulations with N = 100 (we conducted four simulations for each Adisk), giant 
impacts occur after nearly complete sweep up of remnant planetesimals. Therefore, it would 
not be appropriate to apply our analytic estimate to the case of N = 100. Nevertheless, we 
find coincidental good agreements between the averaged e Pjeq 's for simulations with N = 100 
and those from the analytic estimates. 

Finally, let us discuss what will happen if our simulations started with a very small 
size distribution of planete simals. In Figure 8, w e also plot m P)0 ij and m e ff as functions of q 



for the parameters used in IWeidenschilling et al.l (119971 ) (mo = 4.8 x 10 g), who conducted 



planetary accretion simulations in the runaway and oligarchic stages with their multi-zone 
code. For (m Pj0 ii/mo) _9 ~ 2 ^> 1 (this is not the case for our simulations), the minimum 
oligarchic mass is given by 



™ P ,oii = ( ^2vrm -- 2 Sa 2 6 p j j . (36) 

This equation indicates that m P)C ,ii decreases rather strongly with m for small q, as m P)0 ii oc 
m (g+2)/(g+4/3)^ therefore, the minimum oligarchic mass with q = —2.5 for mo = 4.8 x 10 18 g 
is much smaller than ours (m = 3.94 x 10 24 g), and oligarchic bodies started to regulate 
the velocity evolution even when f p is very small (estimated to be ~ 10 -3 ). On the other 
hand, mp >0 u and m e ff for q = —2 are only one order of m agnitude smaller than those for 



1/3 v -l/(g+4/3) 



Run 6. Indeed, in simulations of IWeidenschilling et al.l (119971 ). the maximum mass of remnant 



planetesimals is ~ 10 25 g (whereas oligarchic bodies have masses ~ 10 27 g), suggesting q ~ —2. 
Therefore, from eq. [35] the final equilibrium eccentricities is expected to be smaller than 
those in our simulations only by a factor of 3-4. In the discussion here, we ignored the effect 
of damping due to mutual collisions, which would reduce the eccentricities of planetesimals. 
If the fast clean up of remnant planetesimals happens due to the collisional damping, the 
final planetary system might be unstable, as in the case of a large g. 
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6. Conclusions 

We have conducted direct iV-body simulations of the formation of terrestrial planets 
beginning with a compact planetesimal disk, with the total mass and angular momentum 
being those for the current terrestrial planets. In most of the cases, a planetesimal disk 
results in three planets of similar size, and sometimes a fourth small planet forms around the 
location of Mars. Since a sufficient number of planetesimals remain even after giant impacts 
of protoplanets, orbital eccentricities of the final planets are as small as those of the current 
Earth and Venus. This is a very nice success of our simulations which test the standard 
model for planet formation. 

The final eccentricities of planets are nearly in the equilibrium state for the energy 
partitioning with remnant planetesimals, meaning that the mutual interactions of planets in 
the final state is unimportant in most cases. The final eccentricities of planets depends on the 
initial mass of planetesimals only very weakly, and on the surface density relatively strongly. 
These dependences are interpreted in terms of the effective mass of remnant planetesimals 
(m e ff = (m 2 )/(m)), which determines the strength of the gravitational scattering effects of 
planetesimals. The mass distribution of remnant planetesimals is approximately represented 
by a power law distribution, dn oc m q dm, with q ~ —2, with the upper cut off mass (we 
call it the minimum oligarchic mass) which increases very weakly with the initial mass 
of planetesimals. Therefore, the dependence of the effective mass on the initial mass of 
planetesimals is very weak as well. 

In a few of our simulations, planet-planet interactions are important even at the end of 
simulations and the orbits of final planets can be much more eccentric than for our terrestrial 
planets. This situation seems to appear when the gradient of population of planetesimals 
in the radial direction is large during the giant impact stage. However, the number of our 
simulations is still too small to statistically discuss the conditions that are responsible for 
final planetary orbits. 

We appreciate an anonymous reviewer for useful comments. We are grateful to Derek 
Richardson for providing us with his version of PKDGRAV. We thank Shigeru Ida, Makiko 
Nagasawa, and Eiichiro Kokubo for fruitful discussions. Our simulations have been con- 
ducted with the zBoxl and zBox2 supercomputers at the University of Zurich. We thank 
Doug Potter for his management of the computers. 
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Fig. 1. — Snapshots of Run 6 on the a-e plane. Vertical dashed lines are the inner and outer 
edges of the initial planetesimal disk. The circles represent planetesimals and planets, and 
plotted radius sizes are proportional to the actual radii without artifical enhancement. 
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Fig. 2. — Evolution of the cumulative number vs. mass in units of the initial mass for Run 6. 
The slope of the approximated power-law distribution, dN c oc m q dm, is shown in each panel. 
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Fig. 3. — Evolution of the epicyclic velocity vs. mass for Run 6. The velocity is normalized 
by the Keplerian velocity t>k ep ,mid at a = 0.89AU. The upper and lower horizontal dashed 
lines represent the escape velocities of the largest body and the smallest body, respectively. 
Note that the escape velocity of the smallest body increases by a factor of y/g after 10 5 years 
as we reduce the radius enhancement factor g from 4.3 to 1. 
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Fig. 4. — Time evolution of the averaged quantities for A disk = 0.3AU (Runs 1, lb, 2, and 3). 
From top to bottom, the panels show the number, mass fraction, orbital spacing, and orbital 
eccentricity of planets (m > 2 x 10 26 g), and the ratio of the effective mass of planetesimals 
to the mean mass of planets. In the third panel , the relations between the orbital spacing 
and orbital instability time for e p /h p = 2 and 4 (jYoshinaga et al. I ll999h are shown by lower 
and upper dashed lines, respectively. 
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Fig. 5. — Same as Figure 4 but for the case of A disk = 0.5AU (Runs 5, 5b, 6, and 7). 
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Fig. 6. — Snapshot of all of runs on the a-e plane at 200Myr. The vertical lines represents 
inner and outer edges of the initial planetesimal disk. 
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Fig. 7. — Evolution of the normalized time scales f (upper panel), the velocity of planetes- 
imals v normalized by the escape velocity of planets, and the power-law index, q, for mass 
distribution of planetesimals (lower panel) as functions of mass fraction of planets f v to the 
total mass. T grow , T scat , T damp , and T- in st{g) represent time scales for the growth of planets, 
evolution of the planetesimal velocity due to scattering by planets, damping of eccentricities 
of planets due to the dynamical friction of planetesimals, and orbital instability for multiple 
planet systems, respectively. T grow is slightly shifted downward (by 0.05) to avoid overlapped 
displays of the time scales. 
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Fig. 8. — The minimum oligarchic mass m Pj0 ii (solid lines) and the effective mass of 
planetesimals m c g (dashed lines) verses the power-law index q. The case for our Run 6 
(m = 3.94 x lO ^g, S = 19.1 g cm" 2 , a nd a = 0.89AU), and the case for parameters used 
in simulations of IWeidenschilling et al.l (119971 ) (mo = 4.8 x 10 18 g, S = 16.7 g cm" 2 , and 
a = 1.0AU) are s hown . The latter parameters are also used in IWetherill fc Stewartl (119931 ) 
and llnaba et al.l (1200 ll ) . 
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Fig. 9. — Comparison of the final eccentricities of planets between iV-body simulations (open 
and filled circles; from Table 1 except for iV = 100) and analytic estimates (solid and dashed 
lines; eq. [35]). In analytic estimates, we assume the mass of planets m p to be 0.66M E and 
q = —2 for remnant planetesimals. MVEM and VEM stand for the values for the current 
terrestrial planets with and without Mercury, respectively. 
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Table 1. Initial conditions and final states of simulations 



Run 


Adisk(AU) 


N 


a 




e p (10- 2 ) 


s d (io- 3 ) 


1 


0.3 


1000 


-1 


3 


3.96 


1.55 


lb 


0.3 


1000 


-1 


4 


3.87 


1.51 


2 


0.3 


3000 


-1 


3 


3.87 


1.79 


3 


0.3 


5000 


-1 


3 


3.63 


0.99 


4 


0.3 


3000 


-2 


4 


3.09 


0.99 


■5 


0.5 


1000 


-1 


3 


2.73 


0.95 


5b 


0.5 


1000 


-1 


4 


7.74 


3.49 


G 


0.5 


3000 


-1 


3 


2.98 


0.60 


7 


0.5 


5000 


-1 


3 


2.37 


0.49 


8 


0.5 


3000 


-2 


4 


6.02 


2.97 


MVEM 
VEM 








4 
3 


4.62 
3.41 


1.90 
1.50 



Note. — Parameters Adisk, N, and a stand for the width, number 
of planetesimals, and power-law index for the surface density of ini- 
tial planetesimal disks, respectively, and N p , e p , and Sd stand for the 
number, averaged orbital eccentricity (eq. [3]), and angular momen- 
tum deficit of planets (eq. [8]) at the end of simulations. In Runs lb 
and 5b we switch integrators and reduce g to unity at 5 x 10 4 yr, 
and at 10 5 yr for other runs. MVEM and VEM stand for the current 
terrestrial planets with and without Mercury, respectively. 



